num = 30;
den = conv(conv([1, 0], [0.1, 1]), [0.2, 1]);
G0 = tf(num, den);
margin(G0);

gamma0 = 40; delta = 6;
gamma = gamma0 + delta;
w = 0.01:0.01:1000;
[mag, phase, wout] = bode(G0, w);
n = find(180 + phase - gamma <= 0.1);
wgamma = n(1) / 100;
[mag, phase] = bode(G0, wgamma);
Lhc = -20 * log10(mag);
beta = 10^(Lhc / 20);

w2 = wgamma / 10;
w1 = beta * w2;
numc = [1/w2, 1];
denc = [1/w1, 1];
Gc = tf(numc, denc)
G = Gc * G0
bode(G0, G);
hold on;
margin(G);
beta